bioinformatics.org/bioinformatics.org_software.tar/www.bioinformatics.org/download/ecell2/ecell220setup.exe/{app}/standard/STDR/RapidEquilibriumPReactor.rd

Text File | 2000-03-03 | 5KB | 189 lines

@CLASSNAME: RapidEquilibriumPReactor

@BASECLASS: Reactor
@AUTHOR: Kenta Hashimoto
@EMAIL: [e-mail]
@DATE: 2000 2/2

#%VERSION: E-CELL, Reactor
%VERSION: ecs-1.0, 0.1

@BRIEF_DESCRIPTION: Binding reaction which rapidly reaches equilibrium

@DESCRIPTION: A reactor class for binding reactions which rapidly (within
a step) reaches a state of equilibrium.

This reactor calculates velocity according to the equilibrium constant
inputted by the user. In case where two substrates are already at a
state of equilibrium, velocity is set to zero.

@EQUATION:$$k_{eq} \prod ([S_{k}] - v) =  \prod ([P_{k}] + v)$$

####%SUBSTANCE: max, min

%SUBSTANCE:Substrate, inf, 1
%SUBSTANCE:Product, inf, 1

####%PARAMETER?: NAME,TYPE,DESCRIPTION

%PARAMETER: Keq, Float, $M^{(N_{P}) - (N_{S})}$, Equilibrium constant
@NOTES: ''v'' is velocity. \\ $N_{P}$ is number of Products, $N_{S}$ is number of Substrates.

#######  for .h file ##########

%INCLUDE_FILE_H: StandardHeaders.h

@PRIVATE:
@PROTECTED:
@PUBLIC:

#######  for .C file ##########

#%INCLUDE_FILE_C: RootSystem.C

@INITIALIZE_FUNC:
  int d_Keq = 0;
 // Keq_qua = Keq * N_A^(sum_numPro - sum_numSub) * V^(sum_numPro - sum_numSub)
  for(int i=0; i<numSubstrate(); i++)
    d_Keq -= substrate(i)->coefficient();
  for(int i=0; i<numProduct(); i++)
    d_Keq += product(i)->coefficient();
  Keq *= pow(N_A, d_Keq);
@REACT_FUNC:
{
  Float d = 0;
  int figure = 0;
  Float velocity = 0;
  Float velocity_posi = 1;
  Float velocity_nega = 1;
  int least = 0;
  Float Keq_qua = 1;

  ////
  for(int i=0; i<numSubstrate(); i++)
    for(int j=0; j<substrate(i)->coefficient(); j++)
      Keq_qua *= substrate(i)->substance().supersystem()->volume();
  Keq_qua = Keq / Keq_qua;
  for(int i=0; i<numProduct(); i++)
    for(int j=0; j<product(i)->coefficient(); j++)
      Keq_qua *= product(i)->substance().supersystem()->volume();

  ////
  for(int i=0; i<numSubstrate(); i++)
    for(int j=0 ; j<substrate(i)->coefficient(); j++)
      velocity_posi *= substrate(i)->quantity();
  velocity_posi *= Keq_qua;

  for(int i=0; i<numProduct(); i++)
    for(int j=0; j<product(i)->coefficient(); j++)
      velocity_nega *= product(i)->quantity();

  ////
  if(velocity_posi > velocity_nega){

    for(int i=0; i<numSubstrate(); i++)
      if(substrate(i)->quantity() / substrate(i)->coefficient()
         < substrate(least)->quantity() / substrate(least)->coefficient())
        least = i;

    figure = (int)log10(substrate(least)->quantity()
            / substrate(least)->coefficient());

    d = pow(10L,figure);
    while(figure >= 0){
      velocity += d;

      if(substrate(least)->quantity()
         < velocity * substrate(least)->coefficient()){
        velocity -= d;
        figure--;
    d /= 10;
        continue;
      }

      velocity_posi = 1; velocity_nega = 1;
      for(int i=0; i<numSubstrate(); i++)
    for(int j=0; j<substrate(i)->coefficient(); j++)
      velocity_posi *= substrate(i)->quantity() -
          (velocity * substrate(i)->coefficient());
      velocity_posi *= Keq_qua;
      for(int i=0; i<numProduct(); i++)
    for(int j=0; j<product(i)->coefficient(); j++)
      velocity_nega *= product(i)->quantity() +
        (velocity * product(i)->coefficient());

      if(velocity_posi == velocity_nega)
        break;
      else if(velocity_posi < velocity_nega){
        velocity -= d;
        figure--;
    d /= 10;
      }
      else
        continue;
    }
  }

  else if(velocity_posi < velocity_nega ){

    for(int i=0; i<numProduct(); i++)
      if(product(i)->quantity() / product(i)->coefficient()
         < product(least)->quantity() / product(least)->coefficient())
        least = i;

    figure = (int)log10((double)product(least)->quantity()
                        / product(least)->coefficient());

    d = pow(10L,figure);

    while(figure >= 0){
      velocity -= d;

      if(product(least)->quantity() + (velocity * product(least)->coefficient()) < 0){
        velocity += d;
        figure--;
    d *= 10;
        continue;
      }

      velocity_posi =1; velocity_nega =1;
      for(int i=0; i<numSubstrate(); i++)
    for(int j=0 ; j<substrate(i)->coefficient() ; j++)
      velocity_posi *=
        substrate(i)->quantity() - (velocity * substrate(i)->coefficient());
      velocity_posi *= Keq_qua;
      for(int i=0; i<numProduct(); i++)
    for(int j=0; j<product(i)->coefficient(); j++)
      velocity_nega *=
        product(i)->quantity() + (velocity * product(i)->coefficient());

      if(velocity_posi < velocity_nega){
        continue;
      }
      else if(velocity_posi > velocity_nega){
        velocity += d;
        figure--;
    d /= 10;
      }
      else
        break;
    }
  }
  else
    velocity = 0;

  setActivity(velocity);

  int i = numSubstrate();
  do{
    i--;
    substrate(i)->setQuantity(substrate(i)->quantity() - velocity * substrate(i)->coefficient());
  }while(i!=0);
  i = numProduct();
  do{
    i--;
    product(i)->setQuantity(product(i)->quantity() + velocity * product(i)->coefficient());
  }while(i!=0);
}
@OPTION_C: